{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Example: Obtaining results from the eq and profiles objects\n", "\n", "Here we will take a look at how to access different results stored (or calculated using methods) in the `eq` and `profiles` objects. \n", "\n", "To do this we first need to run a static forward simulation in the MAST-U-like tokamak using previously saved coil currents to generate our results. Of course, the same methods below can be applied after an inverse solve!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Static forward simulation" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "# build machine\n", "from freegsnke import build_machine\n", "tokamak = build_machine.tokamak(\n", " active_coils_path=f\"../machine_configs/MAST-U/MAST-U_like_active_coils.pickle\",\n", " passive_coils_path=f\"../machine_configs/MAST-U/MAST-U_like_passive_coils.pickle\",\n", " limiter_path=f\"../machine_configs/MAST-U/MAST-U_like_limiter.pickle\",\n", " wall_path=f\"../machine_configs/MAST-U/MAST-U_like_wall.pickle\",\n", ")\n", "\n", "# initialise equilibrium object\n", "from freegsnke import equilibrium_update\n", "eq = equilibrium_update.Equilibrium(\n", " tokamak=tokamak,\n", " Rmin=0.1, Rmax=2.0, # radial range\n", " Zmin=-2.2, Zmax=2.2, # vertical range\n", " nx=65, # number of grid points in the radial direction (needs to be of the form (2**n + 1) with n being an integer)\n", " ny=129, # number of grid points in the vertical direction (needs to be of the form (2**n + 1) with n being an integer)\n", " # psi=plasma_psi\n", ") \n", "\n", "# initialise profile object\n", "from freegsnke.jtor_update import ConstrainPaxisIp\n", "profiles = ConstrainPaxisIp(\n", " eq=eq,\n", " paxis=8.1e3,\n", " Ip=6.2e5,\n", " fvac=0.5,\n", " alpha_m=1.8,\n", " alpha_n=1.2\n", ")\n", "\n", "# initialise solver\n", "from freegsnke import GSstaticsolver\n", "GSStaticSolver = GSstaticsolver.NKGSsolver(eq) \n", "\n", "# set coil currents\n", "import pickle\n", "with open('simple_diverted_currents_PaxisIp.pk', 'rb') as f:\n", " current_values = pickle.load(f)\n", "\n", "for key in current_values.keys():\n", " eq.tokamak[key].current = current_values[key]\n", "\n", "# carry out forward solve\n", "GSStaticSolver.solve(eq=eq, \n", " profiles=profiles, \n", " constrain=None, \n", " target_relative_tolerance=1e-9)\n", "\n", "# plot the resulting equilbrium\n", "fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=80)\n", "ax1.grid(True, which='both')\n", "eq.plot(axis=ax1, show=False)\n", "eq.tokamak.plot(axis=ax1, show=False)\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Dashboard of results\n", "\n", "It is worth playing around with the fields/methods in the `eq` and `profiles` objects yourself to see which quantities can be calculated from the plasma equilibrium using built-in functionality. Below, we provide a mini dashboard of different quantities and how to generate them. \n", "\n", "If there are certain quantites that you wish to be added that don't exist within FreeGSNKE at the moment, please do submit a feature request to the Github repository." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### General equillibrium quantites\n", "Here, we display a number of used equilibrium quantites. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plasma quantities\n", "print(rf\"Plasma current: {eq.plasmaCurrent()} [A]\")\n", "print(\"---\")\n", "print(rf\"Poloidal beta (definition 1): {eq.poloidalBeta()}\")\n", "print(rf\"Poloidal beta (definition 2): {eq.poloidalBeta2()}\")\n", "print(rf\"Toroidal beta: {eq.toroidalBeta()}\")\n", "print(rf\"Normalised total beta: {eq.normalised_total_Beta()}\")\n", "print(\"---\")\n", "print(fr\"Plasma internal inductance (li1): {eq.internalInductance1()}\")\n", "print(fr\"Plasma internal inductance (li2): {eq.internalInductance2()}\")\n", "print(fr\"Plasma internal inductance (li3): {eq.internalInductance3()}\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plasma geometry quantities\n", "print(rf\"Minor radius: {eq.minorRadius()} [m]\")\n", "print(\"---\")\n", "print(rf\"Magnetic axis location: {eq.magneticAxis()[0:2]} [m]\")\n", "print(rf\"Geometric axis location: {eq.geometricAxis()} [m]\")\n", "print(rf\"Shafranov shift: {eq.shafranov_shift()} [m]\")\n", "print(rf\"Inner and outer mipdplane (Z=0) radii: {eq.innerOuterSeparatrix(Z=0)} [m]\")\n", "print(rf\"Inner and outer mipdplane (Z=0) radii (alternative): {eq.innerOuterSeparatrix2(Z=0)} [m]\")\n", "\n", "print(\"---\")\n", "print(rf\"LCFS circumference: {eq.separatrix_length()} [m]\")\n", "print(rf\"LCFS area: {eq.separatrix_area()} [m^2]\")\n", "print(rf\"Plasma volume: {eq.plasmaVolume()} [m^3]\")\n", "\n", "print(\"---\")\n", "print(rf\"Aspect ratio: {eq.aspectRatio()}\")\n", "\n", "print(\"---\")\n", "print(rf\"Geometric elongation: {eq.geometricElongation()}\")\n", "print(rf\"Geometric elongation (upper): {eq.geometricElongation_upper()}\")\n", "print(rf\"Geometric elongation (lower): {eq.geometricElongation_lower()}\")\n", "print(rf\"Effective elongation: {eq.effectiveElongation()}\")\n", "\n", "print(\"---\")\n", "print(rf\"Triangularity: {eq.triangularity()}\")\n", "print(rf\"Triangularity (upper): {eq.triangularity_upper()}\")\n", "print(rf\"Triangularity (lower): {eq.triangularity_lower()}\")\n", "\n", "print(\"---\")\n", "s_uo, s_ui, s_lo, s_li = eq.squareness()\n", "print(rf\"Squareness (upper outer): {s_uo}\")\n", "print(rf\"Squareness (upper inner): {s_ui}\")\n", "print(rf\"Squareness (lower outer): {s_lo}\")\n", "print(rf\"Squareness (lower inner): {s_li}\")\n", "\n", "print(\"---\")\n", "L = eq.closest_wall_point()\n", "print(rf\"Point on wall that is closest to the LCFS: {L[0]} [m]\")\n", "print(rf\"Corresponding distance from point to LCFS: {L[1]} [m]\")\n", "\n", "print(\"---\")\n", "print(rf\"Is this a limited plasma?: {eq.flag_limiter}\")\n", "print(rf\"Does the core plasma boundary intersect the wall?: {eq.intersectsWall()}\")\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(rf\"Aspect ratio: {eq.aspectRatio()}\")\n", "print(rf\"Geometric elongation: {eq.geometricElongation()}\")\n", "print(rf\"Effective elongation: {eq.effectiveElongation()}\")\n", "print(rf\"Poloidal beta (method 1): {eq.poloidalBeta()}\")\n", "print(rf\"Poloidal beta (method 2): {eq.poloidalBeta2()}\")\n", "print(rf\"Toroidal beta: {eq.toroidalBeta()}\")\n", "print(rf\"Normalised total plasma beta: {eq.normalised_total_Beta()}\")\n", "print(rf\"Plasma volume: {eq.plasmaVolume()}\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(fr\"Poloidal flux on magnetic axis: {eq.psi_axis} [Webers/2\\pi]\")\n", "print(fr\"Poloidal flux on plasma boundary: {eq.psi_bndry} [Webers/2\\pi]\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# extract the plasma core boundary (specify the number of points you want)\n", "eq.separatrix(ntheta=20) " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# extract the X points (R, Z, psi)\n", "eq.xpt" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# extract the strike points where the separatric intersects the wall, if any (R, Z)\n", "eq.strikepoints()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# here we plot some of the features above on different axes to show where they are located (see code documentation for more details)\n", "\n", "sep = eq.separatrix(ntheta=360)\n", "strikes = eq.strikepoints()\n", "\n", "fig1, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 8), dpi=80)\n", "plt.subplots_adjust(wspace=0.25) # adjust the horizontal space between subplots\n", "\n", "ax1.grid(True, which='both')\n", "eq.tokamak.plot(axis=ax1,show=False)\n", "ax1.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0)\n", "ax1.contour(eq.R, eq.Z, eq.psi(), levels=[eq.psi_bndry], linestyles='solid', colors='r') #, label=\"Separatrix\")\n", "ax1.scatter(strikes[:,0], strikes[:,1], color='k', marker='x', s=40, zorder=2, label=\"Strikepoints\")\n", "ax1.set_aspect('equal')\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "ax1.set_xlabel(r'Major radius, $R$ [m]')\n", "ax1.set_ylabel(r'Height, $Z$ [m]')\n", "ax1.legend(loc=\"upper right\")\n", "\n", "\n", "\n", "\n", "\n", "ax2.grid(True, which='both')\n", "eq.tokamak.plot(axis=ax2,show=False)\n", "ax2.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0)\n", "ax2.scatter(sep[:,0], sep[:,1], color='k', marker='o', s=2, zorder=2, label=\"LCFS (plasma boundary)\")\n", "ax2.scatter(eq.xpt[:,0], eq.xpt[:,1], color='r', marker='x', s=40, zorder=2, label=\"X-points\")\n", "ax2.scatter(eq.magneticAxis()[0], eq.magneticAxis()[1], color='g', marker='x', s=40, zorder=2, label=\"Magnetic axis\")\n", "ax2.scatter(eq.geometricAxis()[0], eq.geometricAxis()[1], color='b', marker='x', s=40, zorder=2, label=\"Geometric axis\")\n", "ax2.set_aspect('equal')\n", "ax2.set_xlim(0.1, 2.15)\n", "ax2.set_ylim(-2.25, 2.25)\n", "ax2.set_xlabel(r'Major radius, $R$ [m]')\n", "ax2.set_ylabel(r'Height, $Z$ [m]')\n", "ax2.legend(loc=\"upper right\")\n", "\n", "\n", "\n", "\n", "\n", "ax3.grid(True, which='both')\n", "eq.tokamak.plot(axis=ax3,show=False)\n", "ax3.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0)\n", "ax3.scatter(sep[:,0], sep[:,1], color='k', marker='o', s=2, zorder=2, label=\"LCFS (plasma boundary)\")\n", "ax3.scatter(eq._sep_Rmin, eq._sep_ZRmin, color='r', marker='x', s=40, zorder=2, label=\"(Rmin, ZRmin)\")\n", "ax3.scatter(eq._sep_Rmax, eq._sep_ZRmax, color='m', marker='x', s=40, zorder=2, label=\"(Rmax, ZRmax)\")\n", "ax3.scatter(eq._sep_RZmin, eq._sep_Zmin, color='g', marker='x', s=40, zorder=2, label=\"(RZmin, Zmin)\")\n", "ax3.scatter(eq._sep_RZmax, eq._sep_Zmax, color='b', marker='x', s=40, zorder=2, label=\"(ZRmax, Zmax)\")\n", "ax3.set_aspect('equal')\n", "ax3.set_xlim(0.1, 2.15)\n", "ax3.set_ylim(-2.25, 2.25)\n", "ax3.set_xlabel(r'Major radius, $R$ [m]')\n", "ax3.set_ylabel(r'Height, $Z$ [m]')\n", "ax3.legend(loc=\"upper right\")\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here we visualise the dr sep quantities.\n", "\n", "This two distances are defined as the radial separation (on the inboard and outboard sides) at the midplane (Z = 0) between flux surfaces passing through the lower and upper X-points." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "dr_seps = eq.dr_sep()\n", "\n", "print(f\"dr_sep_in = {dr_seps[0]} [m].\")\n", "print(f\"dr_sep_out = {dr_seps[1]} [m].\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the plot below, dr_sep_in is the difference between the radial position of the blue and the red flux surfaces at $Z=0$. \n", "\n", "Similarly, dr_sep_out is the difference between the radial position of the red and the blue flux surfaces at $Z=0$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# visualising these quantities\n", "psi_bndry = eq.psi_bndry\n", "xpts = eq.xpt\n", "\n", "# compute absolute difference between each point's ψ and psi_bndry and then\n", "# get indices of the two smallest differences\n", "closest_indices = np.argsort(np.abs(xpts[:, 2] - psi_bndry))[:2]\n", "\n", "# extract the corresponding rows (two X-points closest to psi_bndry, then sort by lowest z coord z-point)\n", "closest_xpts = xpts[closest_indices]\n", "closest_xpts_sorted = closest_xpts[np.argsort(closest_xpts[:, 1])]\n", "\n", "# plot the flux contours for the values of psi_boundary at each X-point\n", "fig1, axes = plt.subplots(1, 3, figsize=(15, 8), dpi=80)\n", "plt.subplots_adjust(wspace=0.25) # adjust the horizontal space between subplots\n", "\n", "axes[0].grid(True, which='both')\n", "eq.tokamak.plot(axis=axes[0],show=False)\n", "axes[0].fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0)\n", "axes[0].contour(eq.R, eq.Z, eq.psi(), levels=[closest_xpts_sorted[0, 2]], linestyles='solid', colors='r')\n", "axes[0].contour(eq.R, eq.Z, eq.psi(), levels=[closest_xpts_sorted[1, 2]], linestyles='dashed', colors='b')\n", "axes[0].scatter(closest_xpts_sorted[0,0], closest_xpts_sorted[0,1], color='r', marker='x', s=40, zorder=2, label=\"Separatrix from lower X-point\")\n", "axes[0].scatter(closest_xpts_sorted[1,0], closest_xpts_sorted[1,1], color='b', marker='x', s=40, zorder=2, label=\"Separatrix from upper X-point\")\n", "axes[0].set_aspect('equal')\n", "axes[0].set_xlim(0.1, 2.15)\n", "axes[0].set_ylim(-2.25, 2.25)\n", "axes[0].set_xlabel(r'Major radius, $R$ [m]')\n", "axes[0].set_ylabel(r'Height, $Z$ [m]')\n", "axes[0].legend(loc=\"upper right\")\n", "\n", "\n", "axes[1].grid(True, which='both')\n", "eq.tokamak.plot(axis=axes[1],show=False)\n", "axes[1].fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0)\n", "axes[1].contour(eq.R, eq.Z, eq.psi(), levels=[closest_xpts_sorted[0, 2]], linestyles='solid', colors='r')\n", "axes[1].contour(eq.R, eq.Z, eq.psi(), levels=[closest_xpts_sorted[1, 2]], linestyles='dashed', colors='b')\n", "axes[1].scatter(closest_xpts_sorted[0,0], closest_xpts_sorted[0,1], color='r', marker='x', s=40, zorder=2, label=\"Separatrix from lower X-point\")\n", "axes[1].scatter(closest_xpts_sorted[1,0], closest_xpts_sorted[1,1], color='b', marker='x', s=40, zorder=2, label=\"Separatrix from upper X-point\")\n", "axes[1].set_aspect('equal')\n", "axes[1].set_xlim(0.18, 0.38)\n", "axes[1].set_ylim(-0.05, 0.05)\n", "axes[1].set_xlabel(r'Major radius, $R$ [m]')\n", "axes[1].set_ylabel(r'Height, $Z$ [m]')\n", "axes[1].set_title(\"Inboard side: dr_sep_in\")\n", "\n", "\n", "axes[2].grid(True, which='both')\n", "eq.tokamak.plot(axis=axes[2],show=False)\n", "axes[2].fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0)\n", "axes[2].contour(eq.R, eq.Z, eq.psi(), levels=[closest_xpts_sorted[0, 2]], linestyles='solid', colors='r')\n", "axes[2].contour(eq.R, eq.Z, eq.psi(), levels=[closest_xpts_sorted[1, 2]], linestyles='dashed', colors='b')\n", "axes[2].scatter(closest_xpts_sorted[0,0], closest_xpts_sorted[0,1], color='r', marker='x', s=40, zorder=2, label=\"Separatrix from lower X-point\")\n", "axes[2].scatter(closest_xpts_sorted[1,0], closest_xpts_sorted[1,1], color='b', marker='x', s=40, zorder=2, label=\"Separatrix from upper X-point\")\n", "axes[2].set_aspect('equal')\n", "axes[2].set_xlim(1.3, 1.5)\n", "axes[2].set_ylim(-0.05, 0.05)\n", "axes[2].set_xlabel(r'Major radius, $R$ [m]')\n", "axes[2].set_ylabel(r'Height, $Z$ [m]')\n", "axes[2].set_title(\"Outboard side: dr_sep_out\")\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Flux quantites\n", "Here, we plot the total poloidal flux $\\psi$ [Webers / $2\\pi$] and its two components $\\psi_p$ (the plasma flux) and $\\psi_c$ (the coil flux): $\\psi = \\psi_p + \\psi_c$. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig1, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 8), dpi=80)\n", "plt.subplots_adjust(wspace=0.5) # adjust the horizontal space between subplots\n", "\n", "ax1.grid(True, which='both')\n", "eq.tokamak.plot(axis=ax1,show=False)\n", "# ax1.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", "ax1.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "im1 = ax1.contour(eq.R, eq.Z, eq.psi(), levels=50) # total psi\n", "ax1.contour(eq.R, eq.Z, eq.psi(), levels=[eq.psi_bndry], colors='r') # psi boundary contour\n", "ax1.set_aspect('equal')\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "ax1.set_title(r\"Total flux, $\\psi$\")\n", "ax1.set_xlabel(r'Major radius, $R$ [m]')\n", "ax1.set_ylabel(r'Height, $Z$ [m]')\n", "cbar = plt.colorbar(im1, ax=ax1, fraction=0.09)\n", "\n", "\n", "ax2.grid(True, which='both')\n", "eq.tokamak.plot(axis=ax2,show=False)\n", "# ax2.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", "ax2.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "im2 = ax2.contour(eq.R, eq.Z, eq.plasma_psi, levels=50) # plasma psi\n", "ax2.set_aspect('equal')\n", "ax2.set_xlim(0.1, 2.15)\n", "ax2.set_ylim(-2.25, 2.25)\n", "ax2.set_title(r\"Plasma flux, $\\psi_p$\")\n", "ax2.set_xlabel(r'Major radius, $R$ [m]')\n", "ax2.set_ylabel(r'Height, $Z$ [m]')\n", "cbar = plt.colorbar(im2, ax=ax2, fraction=0.09)\n", "\n", "\n", "ax3.grid(True, which='both')\n", "eq.tokamak.plot(axis=ax3,show=False)\n", "# ax3.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", "ax3.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "im3 = ax3.contour(eq.R, eq.Z, eq.tokamak_psi, levels=50) # coil psi\n", "ax3.set_aspect('equal')\n", "ax3.set_xlim(0.1, 2.15)\n", "ax3.set_ylim(-2.25, 2.25)\n", "ax3.set_title(r\"Coil flux, $\\psi_c$\")\n", "ax3.set_xlabel(r'Major radius, $R$ [m]')\n", "ax3.set_ylabel(r'Height, $Z$ [m]')\n", "cbar = plt.colorbar(im3, ax=ax3, fraction=0.09)\n", "\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also extract the flux produced by individual coils (or pasive structures) if required. We do this my using the pre-calculated Green's functions, the coil currents, multipliers, and polarities. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "currents = eq.tokamak.getCurrents()\n", "\n", "# calculate the active/passive coil fluxes\n", "\n", "psi_coils = dict()\n", "for i, name in enumerate(currents.keys()):\n", " coil = eq.tokamak.coils_dict[name]\n", " scaling = coil[\"multiplier\"]*coil[\"polarity\"]\n", " greens_matrix = 0.0\n", " \n", " if type(eq._pgreen[name]) is dict:\n", " num_coils = len(eq._pgreen[name])\n", " for i, ind in enumerate(eq._pgreen[name]):\n", " greens_matrix += eq._pgreen[name][ind]*scaling[i*(len(scaling)//num_coils)]\n", " else:\n", " num_coils = 1\n", " greens_matrix = eq._pgreen[name]*scaling[0]\n", "\n", " psi_coils[name] = greens_matrix*currents[name]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here, we plot a few of the fluxes from the coils. You can change the `name` parameters to visualise different coil fluxes." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig1, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 8), dpi=80)\n", "plt.subplots_adjust(wspace=0.5) # adjust the horizontal space between subplots\n", "\n", "name = \"Solenoid\"\n", "ax1.grid(True, which='both')\n", "eq.tokamak.plot(axis=ax1,show=False)\n", "# ax1.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", "ax1.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "im1 = ax1.contour(eq.R, eq.Z, psi_coils[name], levels=50) \n", "ax1.set_aspect('equal')\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "ax1.set_title(rf\"$\\psi_c$ ({name})\")\n", "ax1.set_xlabel(r'Major radius, $R$ [m]')\n", "ax1.set_ylabel(r'Height, $Z$ [m]')\n", "cbar = plt.colorbar(im1, ax=ax1, fraction=0.09)\n", "\n", "name = \"D1\"\n", "ax2.grid(True, which='both')\n", "eq.tokamak.plot(axis=ax2,show=False)\n", "# ax2.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", "ax2.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "im2 = ax2.contour(eq.R, eq.Z, psi_coils[name], levels=50) \n", "ax2.set_aspect('equal')\n", "ax2.set_xlim(0.1, 2.15)\n", "ax2.set_ylim(-2.25, 2.25)\n", "ax2.set_title(rf\"$\\psi_c$ ({name})\")\n", "ax2.set_xlabel(r'Major radius, $R$ [m]')\n", "ax2.set_ylabel(r'Height, $Z$ [m]')\n", "cbar = plt.colorbar(im2, ax=ax2, fraction=0.09)\n", "\n", "name = \"P6\"\n", "ax3.grid(True, which='both')\n", "eq.tokamak.plot(axis=ax3,show=False)\n", "# ax3.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", "ax3.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "im3 = ax3.contour(eq.R, eq.Z, psi_coils[name], levels=50) \n", "ax3.set_aspect('equal')\n", "ax3.set_xlim(0.1, 2.15)\n", "ax3.set_ylim(-2.25, 2.25)\n", "ax3.set_title(rf\"$\\psi_c$ ({name})\")\n", "ax3.set_xlabel(r'Major radius, $R$ [m]')\n", "ax3.set_ylabel(r'Height, $Z$ [m]')\n", "cbar = plt.colorbar(im3, ax=ax3, fraction=0.09)\n", "\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Magnetic fields\n", "\n", "Here, we plot the different magnetic field components over the domain." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig1, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 8), dpi=80)\n", "plt.subplots_adjust(wspace=0.5) # adjust the horizontal space between subplots\n", "\n", "ax1.grid(True, which='both')\n", "eq.tokamak.plot(axis=ax1,show=False)\n", "# ax1.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", "ax1.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "im1 = ax1.contour(eq.R, eq.Z, eq.Br(eq.R, eq.Z), levels=200) \n", "ax1.set_aspect('equal')\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "ax1.set_title(r\"$B_R$ [T]\")\n", "ax1.set_xlabel(r'Major radius, $R$ [m]')\n", "ax1.set_ylabel(r'Height, $Z$ [m]')\n", "cbar = plt.colorbar(im1, ax=ax1, fraction=0.09)\n", "\n", "\n", "ax2.grid(True, which='both')\n", "eq.tokamak.plot(axis=ax2,show=False)\n", "# ax2.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", "ax2.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "im2 = ax2.contour(eq.R, eq.Z, eq.Btor(eq.R, eq.Z), levels=100) \n", "ax2.set_aspect('equal')\n", "ax2.set_xlim(0.1, 2.15)\n", "ax2.set_ylim(-2.25, 2.25)\n", "ax2.set_title(r\"$B_{\\phi}$ [T]\")\n", "ax2.set_xlabel(r'Major radius, $R$ [m]')\n", "ax2.set_ylabel(r'Height, $Z$ [m]')\n", "cbar = plt.colorbar(im2, ax=ax2, fraction=0.09)\n", "\n", "\n", "ax3.grid(True, which='both')\n", "eq.tokamak.plot(axis=ax3,show=False)\n", "# ax3.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", "ax3.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "im3 = ax3.contour(eq.R, eq.Z, eq.Bz(eq.R, eq.Z), levels=100) \n", "ax3.set_aspect('equal')\n", "ax3.set_xlim(0.1, 2.15)\n", "ax3.set_ylim(-2.25, 2.25)\n", "ax3.set_title(r\"$B_Z$ [T]\")\n", "ax3.set_xlabel(r'Major radius, $R$ [m]')\n", "ax3.set_ylabel(r'Height, $Z$ [m]')\n", "cbar = plt.colorbar(im3, ax=ax3, fraction=0.09)\n", "\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# this is just the contribution from the plasma\n", "fig1, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 8), dpi=80)\n", "plt.subplots_adjust(wspace=0.5) # adjust the horizontal space between subplots\n", "\n", "ax1.grid(True, which='both')\n", "eq.tokamak.plot(axis=ax1,show=False)\n", "# ax1.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", "ax1.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "im1 = ax1.contour(eq.R, eq.Z, eq.plasmaBr(eq.R, eq.Z), levels=30) \n", "ax1.set_aspect('equal')\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "ax1.set_title(r\"$B_R$ (from plasma) [T]\")\n", "ax1.set_xlabel(r'Major radius, $R$ [m]')\n", "ax1.set_ylabel(r'Height, $Z$ [m]')\n", "cbar = plt.colorbar(im1, ax=ax1, fraction=0.09)\n", "\n", "\n", "ax2.grid(True, which='both')\n", "eq.tokamak.plot(axis=ax2,show=False)\n", "# ax2.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", "ax2.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "im2 = ax2.contour(eq.R, eq.Z, eq.plasmaBz(eq.R, eq.Z), levels=30) \n", "ax2.set_aspect('equal')\n", "ax2.set_xlim(0.1, 2.15)\n", "ax2.set_ylim(-2.25, 2.25)\n", "ax2.set_title(r\"$B_Z$ (from plasma) [T]\")\n", "ax2.set_xlabel(r'Major radius, $R$ [m]')\n", "ax2.set_ylabel(r'Height, $Z$ [m]')\n", "cbar = plt.colorbar(im2, ax=ax2, fraction=0.09)\n", "\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# this is just the contribution from the active coils\n", "\n", "Br_actives = 0.0\n", "Bz_actives = 0.0\n", "for name in eq.tokamak.coils_list[0:12]:\n", " Br_actives += eq.tokamak[name].Br(eq.R, eq.Z)\n", " Bz_actives += eq.tokamak[name].Bz(eq.R, eq.Z)\n", "\n", "fig1, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 8), dpi=80)\n", "plt.subplots_adjust(wspace=0.5) # adjust the horizontal space between subplots\n", "\n", "ax1.grid(True, which='both')\n", "eq.tokamak.plot(axis=ax1,show=False)\n", "# ax1.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", "ax1.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "im1 = ax1.contour(eq.R, eq.Z, Br_actives, levels=200) \n", "ax1.set_aspect('equal')\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "ax1.set_title(r\"$B_R$ (from coils) [T]\")\n", "ax1.set_xlabel(r'Major radius, $R$ [m]')\n", "ax1.set_ylabel(r'Height, $Z$ [m]')\n", "cbar = plt.colorbar(im1, ax=ax1, fraction=0.09)\n", "\n", "\n", "ax2.grid(True, which='both')\n", "eq.tokamak.plot(axis=ax2,show=False)\n", "# ax2.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", "ax2.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "im2 = ax2.contour(eq.R, eq.Z, Bz_actives, levels=200) \n", "ax2.set_aspect('equal')\n", "ax2.set_xlim(0.1, 2.15)\n", "ax2.set_ylim(-2.25, 2.25)\n", "ax2.set_title(r\"$B_Z$ (from coils) [T]\")\n", "ax2.set_xlabel(r'Major radius, $R$ [m]')\n", "ax2.set_ylabel(r'Height, $Z$ [m]')\n", "cbar = plt.colorbar(im2, ax=ax2, fraction=0.09)\n", "\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# # this is just the contribution from the passive coils (note this is zero because there are no currents in the passives!)\n", "\n", "# Br_passives = 0.0\n", "# Bz_passives = 0.0\n", "# for name in eq.tokamak.coils_list[12:]:\n", "# Br_passives += eq.tokamak[name].Br(eq.R, eq.Z)\n", "# Bz_passives += eq.tokamak[name].Bz(eq.R, eq.Z)\n", " \n", " \n", "# fig1, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 8), dpi=80)\n", "# plt.subplots_adjust(wspace=0.5) # adjust the horizontal space between subplots\n", "\n", "# ax1.grid(True, which='both')\n", "# eq.tokamak.plot(axis=ax1,show=False)\n", "# # ax1.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", "# ax1.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "# im1 = ax1.contour(eq.R, eq.Z, Br_passives, levels=200) \n", "# ax1.set_aspect('equal')\n", "# ax1.set_xlim(0.1, 2.15)\n", "# ax1.set_ylim(-2.25, 2.25)\n", "# ax1.set_title(r\"$B_R$ (from coils) [T]\")\n", "# ax1.set_xlabel(r'Major radius, $R$ [m]')\n", "# ax1.set_ylabel(r'Height, $Z$ [m]')\n", "# cbar = plt.colorbar(im1, ax=ax1, fraction=0.09)\n", "\n", "\n", "# ax2.grid(True, which='both')\n", "# eq.tokamak.plot(axis=ax2,show=False)\n", "# # ax2.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", "# ax2.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "# im2 = ax2.contour(eq.R, eq.Z, Bz_passives, levels=200) \n", "# ax2.set_aspect('equal')\n", "# ax2.set_xlim(0.1, 2.15)\n", "# ax2.set_ylim(-2.25, 2.25)\n", "# ax2.set_title(r\"$B_Z$ (from coils) [T]\")\n", "# ax2.set_xlabel(r'Major radius, $R$ [m]')\n", "# ax2.set_ylabel(r'Height, $Z$ [m]')\n", "# cbar = plt.colorbar(im2, ax=ax2, fraction=0.09)\n", "\n", "\n", "# plt.tight_layout()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# total poloidal magnetic field\n", "\n", "fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=100)\n", "\n", "ax1.grid(True, which='both')\n", "eq.tokamak.plot(axis=ax1,show=False)\n", "# ax1.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", "ax1.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "im1 = ax1.contour(eq.R, eq.Z, eq.Bpol(eq.R, eq.Z), levels=200) \n", "ax1.set_aspect('equal')\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "ax1.set_title(r\"$B_{pol}$ [T]\")\n", "ax1.set_xlabel(r'Major radius, $R$ [m]')\n", "ax1.set_ylabel(r'Height, $Z$ [m]')\n", "cbar = plt.colorbar(im1, ax=ax1, fraction=0.09)\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Coil currents\n", "\n", "Here, we visualise the size of the currents in the poloidal field coils." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# listing the currents may not be that informative \n", "eq.tokamak.getCurrents()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can play around with the colour schemes below to visualise how the current is distributed around the machine. Note here that there is no current in the passive structures but they can be included if present. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from matplotlib.patches import Rectangle, Polygon\n", "from matplotlib.colors import Normalize, TwoSlopeNorm\n", "import matplotlib.cm as cm\n", "\n", "\n", "\n", "# create colormap based on magnitude of currents\n", "currents_array = []\n", "for key in list(eq.tokamak.coils_dict.keys()):\n", " currents_array.append(eq.tokamak[key].current)\n", "\n", "max_curr = np.max(np.abs(currents_array))\n", "# norm = Normalize(vmin=-max_curr, vmax=max_curr) # alternative colorbar\n", "norm = TwoSlopeNorm(vmin=np.min(currents_array), vcenter=0, vmax=np.max(currents_array)) \n", "cmap = cm.bwr\n", "\n", "\n", "# plot\n", "fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=120)\n", "plt.tight_layout()\n", "ax1.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "ax1.grid(alpha=0.5)\n", "ax1.set_aspect('equal')\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "ax1.set_xlabel(r'Major radius, $R$ [m]')\n", "ax1.set_ylabel(r'Height, $Z$ [m]')\n", "\n", "\n", "for name in list(eq.tokamak.coils_dict.keys()):\n", " coil = eq.tokamak.coils_dict[name]\n", " current = eq.tokamak[name].current\n", " color = cmap(norm(current)) # map the current to a color\n", " \n", " # plot active coils (and currents)\n", " if coil[\"active\"]:\n", "\n", " for i in range(0, len(coil[\"coords\"][0,:])):\n", " patch = Rectangle(\n", " (coil[\"coords\"][0,i] - coil[\"dR\"] / 2, coil[\"coords\"][1,i] - coil[\"dZ\"] / 2),\n", " width=coil[\"dR\"],\n", " height=coil[\"dZ\"],\n", " facecolor=color,\n", " edgecolor='k',\n", " linewidth=0.2,\n", " )\n", " ax1.add_patch(patch)\n", " \n", " # plot passive structures (currents are zero here) \n", " else:\n", " patch = Polygon(\n", " coil[\"vertices\"].T, \n", " facecolor=color, \n", " edgecolor=\"k\", \n", " linewidth=0.2\n", " )\n", " ax1.add_patch(patch)\n", "\n", "\n", "# add a colorbar\n", "sm = cm.ScalarMappable(cmap=cmap, norm=norm)\n", "sm.set_array([]) \n", "cbar = fig1.colorbar(sm, ax=ax1, orientation='vertical', fraction=0.09)\n", "cbar.set_label('Current [A]')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Forces on the coils\n", "\n", "We can also extract the radial and vertical forces on the coils." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "eq.printForces()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 1D plasma current density profiles (and others)\n", "Here, we visualise the $p'$ and $FF'$ profiles used in our equilbirium solve. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plot the input p' and FF' profiles\n", "\n", "psi_n = eq.psiN_1D(N=65)\n", "\n", "fig1, (ax1, ax2) = plt.subplots(1, 2, figsize=(15,6), dpi=80)\n", "ax1.grid(zorder=0, alpha=0.75)\n", "ax1.plot(psi_n, profiles.pprime(psi_n), color='k', linewidth=1, marker='x', markersize=2, zorder=10)\n", "ax1.set_xlabel(r'$\\hat{\\psi}$')\n", "ax1.set_ylabel(r\"$p'(\\hat{\\psi})$\")\n", "ax1.ticklabel_format(axis='y', scilimits=(0,0))\n", "\n", "ax2.grid(zorder=0, alpha=0.75)\n", "ax2.plot(psi_n, profiles.ffprime(psi_n), color='k', linewidth=1, marker='x', markersize=2, zorder=10)\n", "ax2.set_xlabel(r'$\\hat{\\psi}$')\n", "ax2.set_ylabel(r\"$FF'(\\hat{\\psi})$\")\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plot q profile\n", "\n", "psi_n = np.linspace(0.01,0.99,65) # values of q at 0 and 1 can be problematic\n", "\n", "fig1, ax1 = plt.subplots(1, 1, figsize=(6,6), dpi=80)\n", "ax1.grid(zorder=0, alpha=0.75)\n", "ax1.plot(psi_n, eq.q(psi_n), color='k', linewidth=1, marker='x', markersize=2, zorder=10)\n", "ax1.set_xlabel(r'$\\hat{\\psi}$')\n", "ax1.set_ylabel(r\"$q(\\hat{\\psi})$\")\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plot fpol\n", "\n", "fig1, ax1 = plt.subplots(1, 1, figsize=(6,6), dpi=80)\n", "ax1.grid(zorder=0, alpha=0.75)\n", "ax1.plot(psi_n, eq.fpol(psi_n), color='k', linewidth=1, marker='x', markersize=2, zorder=10)\n", "ax1.set_xlabel(r'$\\hat{\\psi}$')\n", "ax1.set_ylabel(r\"$fpol(\\hat{\\psi})$\")\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Toroidal current density\n", "\n", "Here, we visualise the toroidal current density inside the core of the plasma. \n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=100)\n", "\n", "ax1.grid(True, which='both')\n", "eq.tokamak.plot(axis=ax1,show=False)\n", "ax1.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "im1 = ax1.contourf(eq.R, eq.Z, eq._profiles.jtor, levels=np.linspace(0.01, np.max( eq._profiles.jtor,), 20)) \n", "ax1.contour(eq.R, eq.Z, eq.psi(), levels=[eq.psi_bndry], colors='r')\n", "ax1.set_aspect('equal')\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "ax1.set_title(r\"Toroidal current density, $J_p$ [A]\")\n", "ax1.set_xlabel(r'Major radius, $R$ [m]')\n", "ax1.set_ylabel(r'Height, $Z$ [m]')\n", "cbar = plt.colorbar(im1, ax=ax1, fraction=0.09)\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# # plot 1D_jtor\n", "# fig1, ax1 = plt.subplots(1, 1, figsize=(6,6), dpi=80)\n", "# ax1.grid(zorder=0, alpha=0.75)\n", "# ax1.plot(eq.psiN_1D(N=101), eq.jtor_1D(N=101), color='k', linewidth=1, marker='x', markersize=2, zorder=10)\n", "# ax1.set_xlabel(r'$\\hat{\\psi}$')\n", "# ax1.set_ylabel(r\"$J_{tor}(\\hat{\\psi})$\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Masking arrays\n", "There are a number of masking arrays that are built during the equilibrium solve that may be useful for plotting or other calculations you may carry out. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig1, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 8), dpi=80)\n", "plt.subplots_adjust(wspace=0.25) # adjust the horizontal space between subplots\n", "\n", "ax1.grid(True, which='both')\n", "# eq.tokamak.plot(axis=ax1,show=False)\n", "# ax1.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor=None, zorder=0)\n", "# ax1.contour(eq.R, eq.Z, eq.psi(), levels=[eq.psi_bndry], linestyles='solid', colors='r') #, label=\"Separatrix\")\n", "ax1.pcolormesh(eq.R, eq.Z, eq.mask_inside_limiter, cmap=\"Greys\")\n", "ax1.set_aspect('equal')\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "ax1.set_xlabel(r'Major radius, $R$ [m]')\n", "ax1.set_ylabel(r'Height, $Z$ [m]')\n", "ax1.set_title(\"Mask of the region allowed to the plasma\")\n", "\n", "\n", "ax2.grid(True, which='both')\n", "# eq.tokamak.plot(axis=ax2,show=False)\n", "# ax2.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor=None, zorder=0)\n", "# ax2.contour(eq.R, eq.Z, eq.psi(), levels=[eq.psi_bndry], linestyles='solid', colors='r') #, label=\"Separatrix\")\n", "ax2.pcolormesh(eq.R, eq.Z, eq._profiles.diverted_core_mask, cmap=\"Greys\")\n", "ax2.set_aspect('equal')\n", "ax2.set_xlim(0.1, 2.15)\n", "ax2.set_ylim(-2.25, 2.25)\n", "ax2.set_xlabel(r'Major radius, $R$ [m]')\n", "ax2.set_ylabel(r'Height, $Z$ [m]')\n", "ax2.set_title(\"Mask within the plasma core boundary\")\n", "\n", "plt.tight_layout()" ] } ], "metadata": { "kernelspec": { "display_name": "freegsnke", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.16" } }, "nbformat": 4, "nbformat_minor": 4 }